function [a_sims_all,chat_sims_all,cstar_sims_all,aphat_sims_all,apstar_sims_all,eta_sims_all,sigma_etasq_sims_all,alpha_star_sims_all,mpc_all,age_sims_all,prior_sims_all]...
    = bc_endsignals_sim_fast(w,param,y_sims_all,eps_eta_mat,reset_shocks,b,mu0_fun, ctilde, stateGrid,stateGrid_step,stateGrid_size,log_c, ex_signals)

% INPUT
    %   param:          parameter vector
    %                   [bta,r,abar,W_cc,kappa,delta,deltax,ybar,sigma_y,
    %                   cstar_ya,cbar,prior_mean_coeff,sigma_c,psi,cov_fun]
    %   y_sims_all:     simulated y                             [sim_periods,nsim]
    %   eps_eta_mat:    eta innovations (epsilon_eta)           [sim_periods,nsim]
    %   reset_shocks:   reset shocks                            [sim_periods,nsim]
    %   b:              borrowing constraint (<=0)              [1,1]
    %   tilde           whether to use ctilde (T) or cstar (F)  boolean

    % OUTPUT
    %   a_sims_all:             simulated values for a              [sim_periods,nsim]
    %   chat_sims_all:          simulated values for chat           [sim_periods,nsim]
    %   cstar_sims_all:         simulated values for cstar          [sim_periods,nsim]
    %   aphat_sims_all:         simulated values for aphat          [sim_periods,nsim]
    %   apstar_sims_all:        simulated values for apstar         [sim_periods,nsim]
    %   eta_sims_all:           simultaed values for eta            [sim_periods,nsim]
    %   sigma_etasq_sims_all:   simulated values for sigma_etasq    [sim_periods,nsim]

    % unpack parameters
    bta = param(1);         
    r = param(2);        
    abar = param(3);
    W_cc = param(4);        
    kappa = param(5);    
    deltax = param(6);
    sigma_c = param(7); 
    psi = param(8); 
    cov_fun = param(9); 
    gam = param(10); 
    sigma_l = param(11);  
    lbar = param(12);  
    delta_eta = param(13);  
    delta_inh= param(14); 
    psi_tau  = 0;     

    [sim_periods, nsim] = size(y_sims_all);

    
    
    a_sims_all          = NaN(sim_periods,nsim);
    chat_sims_all       = NaN(sim_periods,nsim);
    cstar_sims_all      = NaN(sim_periods,nsim);
    aphat_sims_all      = NaN(sim_periods,nsim);
    apstar_sims_all     = NaN(sim_periods,nsim);
    eta_sims_all        = NaN(sim_periods,nsim);
    sigma_etasq_sims_all= NaN(sim_periods,nsim);

    alpha_star_sims_all      = NaN(sim_periods,nsim);
    mpc_all             = NaN(sim_periods,nsim);
    age_sims_all        = NaN(sim_periods,nsim);
    prior_sims_all      = NaN(sim_periods,nsim);

    for j=1:nsim
        %  preallocate
        y_sims           = y_sims_all(:,j);
        a_sims           = NaN(sim_periods,1);
        chat_sims        = NaN(sim_periods,1);
        cstar_sims       = NaN(sim_periods,1);
        aphat_sims       = NaN(sim_periods,1);
        apstar_sims      = NaN(sim_periods,1);
        eta_sims         = NaN(sim_periods,1);
        sigma_etasq_sims = NaN(sim_periods,1);
        sigma_etasq_save = NaN(sim_periods,1);
        eps_eta_sims     = eps_eta_mat(:,j);

        alpha_star_sims = NaN(sim_periods,1);
        age_sims     = NaN(sim_periods,1);
        prior_sims     = NaN(sim_periods,1);
        mpc              =NaN(sim_periods,1);

        %  initial states
        a_sims(1) = abar;
        mu0_ya_hist = [];
        mu0_ya_hist_mpc = [];
        hist_mat = [];
        Linv = [];
        mu_t0_hist = [];


        %  period 1
        if log_c == 0
            [alpha_t, prior_mean ,chat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat, Linv, mu_t0_hist] = bc_update_iid_fast(r,  W_cc, kappa, w*y_sims(1),         a_sims(1),         sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(1),         b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c, psi_tau);
        else
             [alpha_t, prior_mean ,chat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat, Linv, mu_t0_hist] = bc_update_iid_fast(r,  W_cc, kappa, w*y_sims(1),         a_sims(1),         sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(1),         b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c, psi_tau);
            chat = exp(chat); 
        end
        chat_sims(1) = chat;
        aphat_sims(1) = w*y_sims(1)+ (1+r)*a_sims(1)-chat_sims(1);

        eta_sims(1)  = eta_t;
        sigma_etasq_sims(1) = sigma_etasq;
        sigma_etasq_save(1) = sigma_etasq;

        alpha_star_sims(1)=alpha_t;
        prior_sims(1)=prior_mean;
        age_sims(1)=1;


       
        %  period 2+
        for i = 2:sim_periods
            if reset_shocks(i, j) < deltax
                
                
                
                sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1);
                a_sims(i)    = (w*y_sims(i-1) + (1+r)*a_sims(i-1) - chat_sims(i-1))*delta_inh;
                
                %%%% Carry over beliefs properly
            
                sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1)/delta_eta; 
            
                sigma_etasq_hist = sigma_etasq_sims( sigma_etasq_sims < 1e10); 
                eta_hist         = eta_sims( sigma_etasq_sims < 1e10); 
                ya_hist          = y_sims( sigma_etasq_sims < 1e10) + (1+r)*a_sims( sigma_etasq_sims < 1e10); 
                
                [mu_t0_hist, ~] = gp_update_tilde(ya_hist, ya_hist, eta_hist, sigma_etasq_hist, mu0_fun,sigma_c,psi, cov_fun);
            

                mu0_ya_hist = mu0_fun( ya_hist );

            
                hist_mat = [ y_sims( sigma_etasq_sims < 1e10), a_sims( sigma_etasq_sims < 1e10), eta_hist, sigma_etasq_hist];
            
                temp1 = ya_hist - ya_hist';
                Sigma_y1yN = sigma_c^2*exp(-psi*(temp1.^2));
                Sigma22 = Sigma_y1yN + diag(ones(length(ya_hist),1).*(sigma_etasq_hist));
        
                L = chol(Sigma22, 'lower');
            
                Linv = L\eye( length(Sigma22));
                %}
                
                age_sims(i)=1;
                
            else
                %  discount past info: incr variance of eta (decr precision of eta)
                sigma_etasq_sims(1:i-1) = sigma_etasq_sims(1:i-1);
                %  a_t
                a_sims(i)    = w*y_sims(i-1) + (1+r)*a_sims(i-1) - chat_sims(i-1);

                age_sims(i)=age_sims(i-1)+1;

            end

            %  chat_t
            [alpha_t,prior_mean,chat, eta_t, sigma_etasq, mu0_ya_hist, hist_mat,Linv,mu_t0_hist] = bc_update_iid_fast(r, W_cc, kappa, w*y_sims(i), a_sims(i), sigma_c, psi, cov_fun, ex_signals, eps_eta_sims(i),  b, mu0_ya_hist, hist_mat,Linv,mu_t0_hist, ctilde, stateGrid,stateGrid_step,stateGrid_size, log_c, psi_tau);

            if log_c == 1
                chat = exp(chat); 
            end
            
            chat_sims(i)     = chat;
            cstar_sims(i)    = NaN;
            aphat_sims(i)    = w*y_sims(i)+ (1+r)*a_sims(i)-chat_sims(i);
           
        
            eta_sims(i)  = eta_t;  % signal at t
            sigma_etasq_sims(i) = sigma_etasq;  % precision of signal at t
            sigma_etasq_save(i) = sigma_etasq;

            alpha_star_sims(i)=alpha_t;
            prior_sims(i)=prior_mean;

        end

        a_sims_all(:,j)         = a_sims;
        chat_sims_all(:,j)      = chat_sims;
        cstar_sims_all(:,j)     = cstar_sims;
        aphat_sims_all(:,j)     = aphat_sims;
        %apstar_sims_all(:,j)    = apstar_sims;
        eta_sims_all(:,j)       = eta_sims;
        sigma_etasq_sims_all(:,j) = sigma_etasq_save;

        alpha_star_sims_all(:,j)= alpha_star_sims;
        age_sims_all(:,j)= age_sims;
        prior_sims_all(:,j)= prior_sims;
        mpc_all(:,j)=mpc;

    end

end